Revisiting proboscidean phylogeny and evolution through total evidence and palaeogenetic analyses including Notiomastodon ancient DNA

Summary The extinct Gomphotheriidae is the only proboscidean family that colonized South America. The phylogenetic position of the endemic taxa has been through several revisions using morphological comparisons. Morphological studies are enhanced by paleogenetic analyses, a powerful tool to resolve phylogenetic relationships; however, ancient DNA (aDNA) preservation decreases in warmer regions. Despite the poor preservation conditions for aDNA in humid, sub-tropical climates, we recovered ∼3,000 bp of mtDNA of Notiomastodon platensis from the Arroyo del Vizcaíno site, Uruguay. Our calibrated phylogeny places Notiomastodon as a sister taxon to Elephantidae, with a divergence time of ∼13.5 Ma. Additionally, a total evidence analysis combining morphological and paleogenetic data shows that the three most diverse clades within Proboscidea diverged during the early Miocene, coinciding with the formation of a land passage between Africa and Eurasia. Our results are a further step toward aDNA analyses on Pleistocene samples from subtropical regions and provide a framework for proboscidean evolution.


INTRODUCTION
Proboscidea Illiger, 1811 is a mammalian order that includes the extant elephants (Elephas Linnaeus, 1758, and Loxodonta Anonymous, 1827) and a great diversity of extinct species (Shoshani and Tassy, 2005;Gheerbrant and Tassy, 2009). The order originated around 60 million years (Ma) ago in Africa (Gheerbrant, 2009). At present, the oldest proboscidean fossil is Phosphatherium escuilliei Gheerbrant, Sudre & Cappetta, 1996, from Morocco. It comprises cranial and mandibular elements dating to 55 Ma (Gheerbrant, 2009). The evolutionary history of proboscideans is marked by three major radiations. The first occurred during the late Palaeocene/Eocene, with the diversification of primitive proboscideans. The second radiation took place during the early Miocene, with the diversification of ''Gomphotheriidae' ' Hay, 1922, (used sensu lato throughout the text, indicated by brackets) Mammutidae Hay, 1922, andStegodontidae Osborn, 1918, a family within the Elephantoidea. The last radiation took place during the late Miocene/ early Pliocene and resulted in the diversification of Elephantidae Gray (1821), as well part of the superfamily Elephantoidea, including the living elephants (Shoshani and Tassy, 1996;Gheerbrant and Tassy, 2009).
Most of the phylogenetic studies of South American proboscideans are based on morphological data and recovered different possible combinations of evolutionary relationships among these brevirostrine (short jawed), bunodont proboscideans (Shoshani, 1996;Prado and Alberdi, 2008;Ferretti, 2010;Cozzuol et al., 2012;Mothé et al., 2013Mothé et al., , 2016Mothé et al., , 2017a. A molecular study (Buckley et al., 2019) based on proteomic analysis of bone collagen suggested that Notiomastodon would be more closely related to members of the Mammutidae (the family including Mammut americanum, the American true mastodon) than to the Elephantidae, as traditionally found in morphological studies. Also, Nogueira et al. (2021) recently recovered enamel peptides from a N. platensis molar from the Brazilian Intertropical Region, indicating that such peptides could be used as additional data for phylogenetic studies.
In this study, we present almost 3,000 base pairs of mitochondrial ancient DNA (aDNA) obtained from a specimen of N. platensis from Uruguay ( Figure 1) and reconstruct a molecular phylogeny using these iScience Article data together with previously retrieved DNA sequences from seven other members of the Proboscidea. Additionally, we combine molecular and morphological data to obtain a phylogenetic tree spanning all groups of the order Proboscidea.

RESULTS AND DISCUSSION
Ancient DNA authentication and capture efficiency The biomolecular preservation of the sample was clearly near the limit for DNA sequencing. Also, the phylogenetic distance to previously analyzed species of Proboscidea made it challenging to retrieve sufficient sequence data for analysis. Despite extensive efforts, we were only able to recover slightly less than 3 kb of mitochondrial sequence. The poor preservation of the recovered DNA is exemplified by the characteristics of the reads identified as deriving from Notiomastodon. First, the read length is very short, with a mean of 34 bp and less than 65 reads longer than 50 bp ( Figure S1B). Second, despite treatment with uracil glycosylase, the ends of the fragments show deamination rates of about 30% ( Figure S1C). We do not see any evidence for the presence of nuclear mitochondrial DNA (NUMTs; e.g., heterozygosity) and could exclude the presence of frameshifts or stop codons in open reading frames.
In addition to poor DNA preservation, previous work has shown that hybridization capture fails to recover sequences that are too divergent from the bait sequences (Paijmans et al., 2016(Paijmans et al., , 2020. Since no sequences from close relatives are available, we used baits designed from members of the Elephantidae as well as the American mastodon. As is expected under such circumstances, we preferentially captured relatively conserved regions ( Figure 2). Interestingly, for most of these regions, coverage was relatively high (Figure S1A). This is a direct result of the multiple libraries that were used to build the consensus sequence. Each library contains unique molecules and is therefore deduplicated separately, resulting in stacked reads in regions covered by more than one library. This suggests that it was more the genetic distance between baits and target, rather than a lack of mitochondrial target DNA in the extract, that prevented retrieval of a larger fraction of the mitochondrial genome. It has been shown that capturing mitochondrial genomes using a divergent bait sequence can lead to a loss of regions with more sequence divergence in favor of those with less sequence divergence (Paijmans et al. 2016(Paijmans et al. , 2020Portik et al., 2016, Peñ alba et al., 2014. Additionally, variation in coverage has been observed in the capture of divergent species in previous studies (Paijmans et al., 2020). Divergence between bait and target does not seem to be the sole driver of this effect. Although the sequence divergence is expected to be significant between Notiomastodon and other proboscideans (divergence between the American mastodon and the woolly mammoth is between 3.5% and over 15%; Figure 2), successful capture has been reported for more divergent sequences (40%) in previous papers (e.g., Li et al., 2013). A previously proposed explanation for this large variation in coverage has been the bias introduced by amplification cycles or by the double-capture strategy (Paijmans et al., 2020).
Shotgun libraries contained from zero detectable mitochondrial reads to 0.0002% (Table S1). Despite the phylogenetic distance, the capture strategy was partially successful, leading to an at least 3.53 increase in mitochondrial content (endogenous mitochondrial DNA content for enriched samples used in constructing the consensus sequence was between 0.0009% and 0.0043%; both times mapped to Loxodonta africana; Table S1). It is unclear if changes in capture design would have resulted in higher enrichment, because highly divergent sequences have so far not been shown to profit from, e.g., lowering of hybridisation temperature (Paijmans et al., 2016). To the contrary, lowering the hybridisation temperature might prove counterproductive for libraries containing high levels of contamination (Paijmans et al., 2016). Additionally, due to the high levels of exogenous DNA in our libraries, both the relaxed and the iterative mapping approach we tried failed to yield interpretable results. As to our knowledge, no DNA sequences are so far available for the ''Gomphotheriidae'', our analyses, nevertheless, represent another major step toward extending proboscidean paleogenetics from cold and temperate regions into warmer regions, where a number of recently extinct species and families are awaiting molecular analyses. There have been several recent publications reporting aDNA sequence data from tropical regions (e.g., Kehlmaier et al., 2017Kehlmaier et al., , 2021Woods et al., 2018), but all these investigated samples are only a few thousand years old. While the geographic area in Uruguay itself is defined as humid sub-tropical, the sample age of 35 ka greatly extends the time range of previous studies. The fossils' depth-of-burial and the anoxic geochemical conditions most likely prevented total DNA loss that would be expected in warm climates.
Moreover, the high coverage of part of the mitochondrial genome suggests that sufficient ancient DNA is potentially available in samples from this site, at least for reconstructing mitochondrial genomes. As it may ll OPEN ACCESS iScience 25, 103559, January 21, 2022 3 iScience Article be difficult to design more efficient capture approaches, future studies might either target specific skeletal elements that are expected to have a better DNA preservation (e.g. Alberti et al., 2018, Pinhasi et al., 2015 or try to improve the ratio of endogenous to exogenous DNA during extraction and library prepration (e.g. Essel et al., 2021). As the approach using bleach prior to extraction (Korlevi c et al., 2015) was not successful in our case (see Table S1, libraries STE_B1_cap2 and STE_B2), targeting specific skeletal elements is more likely to yield improved DNA recovery.

Phylogenetic reconstruction and divergence time estimates
The 2,929-bp alignment used for calculating the calibrated tree still contained 152 polymorphic positions, 54 of which are parsimony informative. Our analyses place Notiomastodon as a sister taxon to the Elephantidae, a relation also found by Mothé et al. (2016) using morphological data for a ''Gomphotheriidae'' review. Despite the comparatively short length of the DNA sequence used for phylogenetic analyses iScience Article (Rohland et al., 2007), support for a sister group relationship of Notiomastodon and Elephantidae is strong with a posterior probability of 1 ( Figure 3) and a bootstrap value of 100 ( Figure S1D). Moreover, the topologies recovered both for the seven Elephantidae species included in the Bayesian analysis and the 35 taxa in the maximum likelihood approach are identical to those obtained with full mitochondrial genomes (Figure S1D), except for the placement of Clade 1 mammoths within mammoth diversity, supporting the reliability of the topologies recovered with our dataset.
The phylogenetic position we recover-Notiomastodon being more closely related to Elephantidaeagrees with several morphological studies that recovered this South American proboscidean within the clade of trilophodont gomphotheres (Gheerbrant and Tassy, 2009;Shoshani, 1996;Shoshani et al., 1998Shoshani et al., , 2007 but disagrees with a recently published paleoproteomic study on the basis of collagen sequences (Buckley et al., 2019). In the latter study, Notiomastodon is placed closer to the Mammutidae, an unprecedented phylogenetic hypothesis. However, collagen is known to be a highly conserved biomolecule and provides only limited taxonomic resolution (Cappellini et al., 2018), an observation the authors emphasize as potentially limiting their analysis (Buckley et al., 2019).
Our molecular dating allows an estimation of the divergence time between Notiomastodon and the obtained sister lineage, the Elephantidae. We recover with about 13.5 Ma a much deeper divergence for Notiomastodon versus Elephantidae than for the most recent common ancestor of the Elephantidae, which is estimated to be around 8 Ma, again in agreement with previous results (Gheerbrant and Tassy, 2009;Mothé et al., 2016; Table S2). We tested the effect of different combinations of all three, two, or one internal fossil calibrations as well as the effect of the choice of calibration schemes and prior distribution. Exclusion of the oldest calibration point (split between the American mastodon and all other analyzed genera; Elephantid-Mastodon in Table S2) shifts the estimated divergence between Notiomastodon and the Elephantidae by 2-2.5 Ma to 11-12 Ma, which is to be expected, because this is the only calibration point that is older than the split between Notiomastodon and the Elephantidae. When all three calibration points were used, differences in calibration scheme (broad vs. narrow ages for fossil calibrations) as well as the choice of prior distribution (uniform, log-normal, normal) are not pronounced, with normal prior distributions being the most consistent between the two calibrations schemes (Table S2), whereas both log-normal and uniform prior calibrations differ by 300-400 ka between schemes. There was no apparent directional effect visible  Table S2 for other combinations of fossil calibrations and prior distributions). Blue bars represent 95% highest posterior probability intervals for node ages. Node support is given as Bayesian posterior probability ll OPEN ACCESS iScience 25, 103559, January 21, 2022 5 iScience Article in the choice of prior distribution, e.g., using a log-normal prior distribution made the split younger when using all three as well as only the oldest calibration point, but older when using either the split between the African and the Eurasian lineages (Loxodonta-Eurasian in Table S2) or between the Asian elephants and the mammoths (Asian-mammoth in Table S2).
Combining molecular and morphological data in a total evidence analysis For the majority of fossil proboscideans, no aDNA sequence data are available because most species are beyond the age of the oldest aDNA presently recovered, an approximately 1.5-million-year-old mammoth specimen (van der Valk et al., 2021). Although there has been enormous progress analyzing ancient proteome data, reaching greater time depth than aDNA (Wadsworth and Buckley, 2014;Demarchi et al., 2016), data availability is limited and many extinct proboscidean species also exceed the time frame of proteome analyses, which is currently 2 million years for fossil enamel from temperate or tropical regions (Cappellini et al., 2019). Therefore, to obtain a better understanding of the relationships within the Proboscidea, we performed a total evidence analysis, combining the available DNA sequence data with a large morphological character matrix. The results of this analysis show a topology similar to previously published phylogenetic hypotheses of proboscidean evolution (Figures 4 and S2;Gheerbrant and Tassy, 2009;Shoshani, 1996;Shoshani et al., 2007Shoshani et al., , 1998). The recovered phylogeny shows an early diverging lineage related to the genus Mammut (family Mammutidae, the true mastodons) and the existence of three large groups concordant with the Amebelodontidae, ''Gomphotheriidae'', and Elephantoidea. The results indicate relatively low support for the divergence times of these three commonly recognized groups, which do not allow to confidently resolve the relationships among them, with ''Gomphotheriidae'' still being paraphyletic, whereas the other two represent monophyletic groups. However, we performed additional analyses excluding two taxa that can be considered conflicting due to their probably paraphyletic nature-Gomphotherium and Choerolophodon (see Shoshani, 1996;Tassy, 1996;and Li et al., 2019 for discussions on the subject). The results obtained with this modified dataset showed an identical topology for the remaining taxa and much more support for the basal nodes, indicating that the topology obtained in the main test can be tentatively considered for further analyses ( Figure S3). Regarding Elephantoidea, the analysis recovered Stegodon within the Elephantinae, although node support is low. This result is in contradiction with morphological parsimony results (Shoshani, 1996) that show Stegodon as forming a clade with Stegolophodon in the family Stegodontidae. In our results, Stegolophodon is placed as the most basal taxon in Elephantoidea, also contradicting parsimony results where it was placed within the clade. Finally, in relation to the clade including the trilophodont gomphotheres from Asia and the New World, the results are largely consistent with previous parsimony results (Shoshani, 1996), showing a clade comprising Gnathabelodon, Eubelodon, Rhynchotherium, Sinomastodon, Cuvieronius, Stegomastodon, and Notiomastodon. Our results show Notiomastodon as closely related to Stegomastodon, contrary to recent results by Mothé et al. (2016), but in agreement with the results of Shoshani (1996). However, it must be mentioned that the morphological data employed for the analysis are largely based on Shoshani (1996); new molecular and morphological data will be needed to better resolve the clade's phylogenetic relationships. Nevertheless, despite the long and complicated taxonomic history of these proboscideans, the North American species of Stegomastodon and the South American Notiomastodon are both valid and distinct taxa, according to recent morphological reviews (Mothé et al., 2012(Mothé et al., , 2017a(Mothé et al., , 2017bMothé and Avilla, 2015).
Additionally, we used this combined dataset for dating the proboscidean phylogeny. Regarding the obtained divergence times, the results were largely consistent with the molecular-only analysis. Interestingly, the results showed that the divergences of the three major groups, Amebelodontidae, ''Gomphotheriidae'', and Elephantoidea, occurred during the early Miocene, ca. 18.9 Ma and 17.1 Ma (with a maximum span from 21.2 Ma to 15.1 Ma considering the 95% confidence range of both events). The later divergence (''Gomphotheriidae''-Elephantoidea) date is slightly older than the value obtained in the molecular-only analysis (13.5 Ma;

Diversification and historical biogeography analysis
The diversification analysis conducted in RPANDA identified a single significant shift in the diversification dynamics located in the branch leading to the Elephantida (''Gomphotheriidae'' and Elephantoidea) node ll OPEN ACCESS 6 iScience 25, 103559, January 21, 2022 iScience Article and immediately before the divergence of the three major groups (Figure 3; BIC random /BIC test = 4.36). Interestingly, this period is widely acknowledged as the time that a connection between Africa and Eurasia formed (16-19 Ma;Rö gl, 1998;Harzhauser et al., 2007). This land passage, commonly known as ''The Gomphotherium Landbridge,'' allowed the passage of many mammals that were previously confined to Africa into Eurasia, including the Proboscidea (Sen, 2013;Rage and Gheerbrant, 2020). Regarding the New World brevirostrine gomphotheres, the node representing the last common ancestor of Cuvieronius, Stegomastodon, and Notiomastodon was dated to 5.49 Ma, similar to previous time estimates based on species' stratigraphic ranges (Mothé et al., 2016).
The historical biogeography analysis (Figures 4, S4, and S5) suggests that proboscideans may have left Africa only three times in phylogenetically distant clades (Deinotherium or Deinotheriidae, Mammutidae, and Elephantida). After the connection of Africa to Eurasia, a rapid global dispersal event is observed, which is related to the diversification shift obtained in the RPANDA analysis. These results match the fossil record for Proboscidea, where the first undisputed members (Prodeinotherium and Gomphotherium) of the clade outside Africa are recorded in Pakistan during the Late Oligocene-Early Miocene (Sen, 2013). iScience Article Dispersal from Eurasia back into Africa happened even less often, potentially in the lineage ancestral to Elephantidae. There is also some evidence for additional dispersal between Africa and Eurasia leading to the hybrid origin of the European straight-tusked elephant (Palkopoulou et al., 2018), but when and in which direction this occurred is unknown. At least two large vicariance events are observed, one among the ancestors of the extant genera Elephas and Loxodonta and another marking the arrival of the ancestor of Notiomastodon to North America. Proboscideans apparently reached the Americas independently at least three times, with the ancestors of the American mastodon (Mammutidae), the shoveltusked Amebelodontidae, and the mammoths (Elephantidae). There is also evidence for repeated Pleistocene mammoth dispersal between Eurasia and North America (Debruyne et al., 2008;van der Valk et al., 2021).
Finally, the dispersal to South America was modeled differently when considering the two age constraints for the formation of the Isthmus of Panama. When the age constraint was set to 3.5 Ma, two dispersal events were recorded, one for Cuvieronius and the other for Notiomastodon ( Figure S4). When the constraint was set to 9 Ma, a single dispersal to South America as early as 5.5 Ma was obtained for the common ancestor of Cuvieronius, Notiomastodon, and Stegomastodon ( Figure S5). The earliest undisputed record of a proboscidean in South America is at 2.5 Ma in Argentina (Reguero et al., 2007). This age is consistent with using the late age constraint. Nonetheless, the controversial Amahuacatherium from the Late Miocene of Peru could represent an early dispersal event coinciding with an early ''Great American Biotic Interchange'' (GABI) that is compatible with the 9-Ma constraint (Campbell et al., 2000(Campbell et al., , 2009. However, several arguments regarding the validity of its diagnostic features were recently published (Alberdi et al., 2004;Ferretti, 2008;Mothé and Avilla, 2015). Consequently, the validity of Amahuacatherium is still questionable, with some authors considering the specimen as possibly belonging to Notiomastodon (Mothé and Avilla, 2015) and probably originating from the Pleistocene (Alberdi et al., 2004;Ferretti, 2008;Lucas, 2013).
Considering this, a late arrival scenario is more plausible for South American gomphotheres. The GABI was probably a complex, multi-phase process, with many intervening factors of differential importance in each clade (O'Dea et al., 2016;Jaramillo, 2018). In the case of proboscideans, southward emigration must have occurred at a late stage. As long-distance swimming is known for modern elephants (Johnson, 1980) and presumed for other taxa within the clade, gomphothere's late arrival in South America may be due more to unsuitable environmental conditions before the Pliocene and Pleistocene (Jaramillo, 2018) than to lack of a continuous land bridge.

CONCLUSIONS
Our analysis of aDNA sequences from a 35,000-year-old fossil from a Uruguayan site marks a further step forward in the retrieval of Pleistocene aDNA from humid subtropical regions. Given the numerous and often unique, extinct species from warm climatic regions, our results portend new possibilities for a better understanding of these species' evolutionary history. Our results also show that because living and late Pleistocene representatives (i.e., those that are within the temporal window of DNA preservation) are a relatively small proportion of the taxa of any animal group, total evidence analysis is a useful and underappreciated tool for revealing evolutionary processes and events on a larger temporal scale. Our application of this approach to the order Proboscidea allows insights into their evolution dating back almost 40 Ma and resolves most of their divergence events. However, some internal phylogenetic relationships still remain unresolved in our analyses, especially regarding several representatives of the ''Gomphotheriidae''. Nevertheless, our resolution of the evolutionary and biogeographical history of N. platensis, the unique, endemic proboscidean species of South America, substantially advances our understanding of the natural history of this group in the New World.

Limitations of the study
Our study is limited by representing a single individual of N. platensis and by the fact that we were only able to recover partial mitochondrial data. Additionally, phylogenies recovered from mtDNA can differ from the species tree. The divergence times we recover using mitochondrial DNA are estimates, and new fossil finds could alter our calibration scheme and therefore change our current divergence time estimates. Additionally, a limitation of our study in regard to the total evidence analysis is the low diversity of bunodont proboscideans with recovered aDNA-only Notiomastodon-considering their high diversity and the wide geographic distribution of the group. Considering the approach to the historical biogeography analysis

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
The authors declare no competing interests. Approximately 1 m of alluvium fills the 11-meter long (W-E) basin. The basal stratum, Bed 1, lies unconformably on the Cretaceous sandstone and is green muddy sand that is not fossiliferous. Lying unconformably above is the highly fossiliferous Bed 2, which consists of basal, muddy sandy gravels that fine-upward into muddy sands. The uppermost stratum, Bed 3, forms the modern stream bed and comprises modern, reworked Bed 2 sediments and occasional Bed 2 bones. The age estimate for Bed 2 fossils was reassessed using improved collagen purification techniques, i.e., ultrafiltration, which isolates >30,000 molecular weight (>30 kDa) gelatinized collagen ( Figure S6). These new 14 C dates on bone from Lestodon, Smilodon, and Notiomastodon provide a more accurate date range of 29.9 to 30.7 ka RC years for the bone bed and a direct age of 30,510 G 240 RC yr (UCIAMS-142845) for the Notiomastodon platensis tooth used for aDNA analyses ( Figure S6). Until a larger number of taxa are dated, the calibrated age range of Arroyo del Vizcaíno fossils is assumed to be 34.8 to 35.6 ka cal BP.

Specimen description
The Notiomastodon tooth analysed in this report was a single, isolated 1 st right upper molar (RM 1 ) from the middle of Bed 2, completely surrounded by a fine matrix and in proximity to a sloth tibia. The tooth is one of three Notiomastodon skeletal elements found in Bed 2.
The tooth was very use-worn ( Figure 1B; wear stage 4 of Simpson and Paula-Couto, 1957;Mothé et al., 2010). It presents the typical dental features of South American proboscideans, such as the bunodont morphology (blunt and rounded pairs of cusps, the lophs), three lophs (since it is a molar from intermediate dentition), and a distal cingulum composed of several small accessory cusps ( Figure 1B). The roots were not preserved. The crown structures, e.g., the main and secondary cusps and the accessory cusps have largely eroded during mastication (Shoshani and Tassy, 1996;Mothé and Avilla, 2015). Enamel traces are still present on the labial and lingual borders, and on the labial interloph regions. According to the type of molar (RM 1 ) and its wear stage, we inferred that it belonged to an adult individual of Notiomastodon platensis, possibly around 20 years of age (Mothé et al., 2010). The specimen is archived at the SAUCE-P Laboratory and Collection under the collection number CAV 499, Sauce, Uruguay.

DNA extraction and library preparation
All palaeogenetic work was carried out in designated ancient DNA facilities at the University of Potsdam. Pieces of the dentine were ground to powder using mortar and pestle or a microdismembrator at a frequency of 30 Hertz for 10 sec. Nineteen extracts were prepared, each using approximately 50 mg total ll OPEN ACCESS iScience 25, 103559, January 21, 2022 15 iScience Article of dentine powder, and extracted following the protocol of Dabney et al. (2013). Six of these extracts were obtained after the powder had been treated with 1 mL of 1% sodium hypochlorite for 15 mins prior to extraction (Korlevi c et al., 2015). All samples were incubated overnight in 1 mL extraction buffer (0.45 M EDTA, 0.25 mg/mL Proteinase K) at 37 C under constant rotation. Next, undigested material was pelleted through centrifugation and the supernatant was transferred into 13 mL of binding buffer (5 M guanidine hydrochloride, 40% isopropanol, 0.05% Tween-20, and 90 mM sodium acetate), which was then passed through QIAGEN MinElute columns fitted with a reservoir (Zymo-Spin V). The mix of supernatant and binding buffer for samples STE_CO3 and STE_CO4, as well as STE_UV5 and STE_UV6 were run through the same purification column, respectively, in an attempt to get higher concentrations in the final, purified extracts (Table S1). This was followed by two subsequent wash steps with PE buffer (QIAGEN) and a dry spin of 1 min at 13,000 rpm. The purified DNA was eluted in 25 mL TET buffer (10 mM Tris-HCl, 1 mM EDTA, 0.05% Tween-20).
All seventeen extracts were converted into single-stranded libraries including UDG treatment to reduce DNA damage (Gansauge and Meyer, 2013). To remove residual phosphate groups, one Unit of FastAP was used. Then the DNA was denatured at 95 C for 2 min, resulting in single-stranded DNA. Adaptor CL78 was now ligated to the 3 0 end in an 80 mL reaction containing 20% (vol/vol) PEG-4000, 0.125 mM CL78, and 2.5 U/mL Circligase II, which was incubated overnight at 60 C. The DNA was immobilised on streptavidin covered magnetic beads (Dynabeads MyOne C1) and extension primer CL9 was annealed to the complementary CL78 adaptor. To fill in the second strand, a 50 mL reaction with the following reagents was used: 13 isothermal amplification buffer, 250 mM of each dNTP, 2 mM CL9 extension primer, and 0.48 U/mL Bst 2.0 polymerase. This was followed by a 100 mL reaction to achieve blunt-ended molecules that contained 13 Buffer Tango, 0.025% (vol/vol) Tween 20, 100 mM of each dNTP, and 0.05 U/mL T4 DNA polymerase. In a 100 mL reaction containing 13 T4 DNA ligase buffer, 5% (vol/vol) PEG-4000, 0.025% (vol/ vol) Tween 20, 2 mM double-stranded adaptor, and 0.1 U/mL T4 DNA ligase, the second, double-stranded adaptor (CL53/CL73) was ligated to the molecules. As before, we incubated the mixture at 95 C for 1 min to denature the DNA molecule and the strand complementary to the original single-stranded molecule was eluted in 25 mL of TET buffer. Libraries were amplified and double-indexed in 80 mL reactions containing 13 AccuPrime Pfx reaction mix, 10 mM each of P5 and P7 indexing primers, and 0.025 U/mL AccuPrime Pfx polymerase. The optimal number of cycles was determined by qPCR prior to amplification. This was done in 10 mL reactions containing 13 SYBR green qPCR master mix, 0.2 mM each of IS7 and IS8 amplification primers, and 1 mL of a 1:20 dilution of the unamplified library.

Enrichment and sequencing
Twelve libraries were arbitrarily chosen for enrichment of mitochondrial DNA (mtDNA) using the same capture probes as described in Meyer et al. (2017), but following the protocol described in Gonzalez-Fortes and Paijmans (2019). Library and bait were pooled in a 10:1 ratio, mixed with 2 mM of each blocking oligonucleotide, 13 HI-RPM hybridization buffer and 13 Blocking Agent. We performed two subsequent rounds of enrichment, each time starting at 95 C for 5 min and then cooling down to 65 C at 0.1 C/sec, holding the temperature at 65 C for 24 hours. The capture reactions included extraction and library blanks as well as additional capture blanks.
The 12 enriched libraries and 5 shotgun libraries were sequenced on an Illumina NextSeq 500 as described in Paijmans et al. (2017), either in single end mode (SE, 1 x 75 cycles) or in paired end mode (PE, 2 x 75 cycles). To check for the presence of contamination, extraction and library blanks were included both in shotgun sequencing as well as after enrichment.

Data processing and consensus calling
Adapter sequences and low-quality bases (-q 30) were trimmed using the program cutadapt 1.10 (Martin, 2011). Sequences shorter than 25 bp were removed (-m 25). FLASH 1.2.11 (Mago c and Salzberg, 2011) was used to merge paired end sequencing data with a minimum overlap of 15 bp (-m 15) and a maximum of 10% difference allowed (-x 0.1).
Reads were mapped against the genome of the African Savannah elephant Loxodonta africana (loxAfr3, GenBank Assembly ID: GCA_000001905.1), the mitochondrial genome of the same species (GenBank: NC_000943) as well as mitogenomes of Mammuthus primigenius (GenBank: NC_007596) and Mammut ll OPEN ACCESS

Competitive mapping to human mtDNA
To exclude the presence of human derived reads, we remapped all reads previously mapped to the woolly mammoth reference sequence (GenBank: NC_007596) to a concatenated sequence of the former and the human mitochondrial reference sequence (GenBank: NC_012920) in a competitive mapping approach (Feuerborn et al., 2020), which allows for the identification of human DNA in our mapped reads. 74 reads (3.35%) mapped better to the human reference sequence, but we could exclude any human-derived substitutions in our Notiomastodon consensus sequence due to our strict consensus calling approach.
We used the same alignment as above, added our Notiomastodon consensus sequence and removed all columns with missing data resulting in an alignment of 2940 bp. We then calculated pairwise divergence in non-overlapping 5bp windows and again plotted the data as described above ( Figure 2B).

Molecular phylogenetic analysis
The consensus sequence for Notiomastodon platensis was aligned to seven proboscidean mitochondrial reference sequences as described above and columns with missing data were removed from the alignment, resulting in a final length of 2929 bp.
A Bayesian analysis was performed in BEAST 1. 8.4 (Drummond et al., 2012), in order to estimate the divergence time between Notiomastodon platensis and the Elephantidae. An appropriate partitioning scheme from all possible combinations of available genes, two rRNAs, tRNAs, and part of the D-loop was selected ll OPEN ACCESS iScience 25, 103559, January 21, 2022 iScience Article using PartitionFinder 1.1.1 under the Bayesian Information Criterion (Lanfear et al., 2012). Because of the incomplete coverage of our sequence, not all genes and tRNAs were available and we did not distinguish individual coding positions for genes. This resulted in only one partition with the TrN+I substitution model. We applied a speciation birth-death tree prior and tested for rate variation by using a lognormal relaxed clock. Since the variation of the substitution rate among lineages was found to abut zero, we reran the analysis with a strict clock and applied this to all successive runs.
To calibrate the molecular clock, we ran several analyses using either narrow fossil calibrations (Brandt et al., 2012) or broad fossil calibrations (Rohland et al., 2010). Detailed justifications for the broad fossil calibration scheme can be found in Supplementary Text S5 in Rohland et al. (2010) and are mainly based on Sanders et al. (2010). Narrow fossil calibrations are based on the broad fossil calibration scheme but exclude questionable taxa and assume monophyly for Elephas and Loxodonta, which is not given within the broad fossil calibration scheme (Rohland et al., 2010;Brandt et al., 2012). We calibrated up to three internal nodes, using all combinations of one, two or all three calibration points. Prior distributions for fossil priors were drawn from uniform, lognormal or normal distributions (see Table S2). Additionally, we used a maximum-likelihood approach on a bigger data set to reconstruct the phylogenetic position of Notiomastodon. Thirty-five proboscidean mitochondrial sequences (see Figure S1D for GenBank accession numbers) were aligned as described before. From this alignment (excluding Notiomastodon) we calculated a maximum-likelihood tree using RAxML-HPC 8.2.4 (Stamatakis, 2014) CIPRES black box version on the CIPRES Science Gateway (Miller et al., 2010), with the default GTR+CAT substitution model in place. This phylogeny of complete mitochondrial genomes (16,026 bp) was used in comparison to a phylogeny that included the incomplete sequence of Notiomastodon and therefore reduced the length of the alignment. We added the sequence of Notiomastodon platensis to the alignment using MAFFT as described before and removed missing data. This resulted in an alignment of 2548 bp length, which is shorter than the original consensus sequence of Notiomastodon. This reduction in the length of the alignment is caused by the presence of missing data in many published woolly mammoth mitochondrial sequences. A maximum-likelihood tree was calculated as described before and the recovered topologies were compared for consistency.
Total evidence phylogenetic analysis. We performed a total evidence approach by combining the new molecular data with published morphological data for a wide range of proboscidean fossil taxa. We updated the morphological matrix of Shoshani (1996) and Shoshani et al. (2006), and added recent data for several of the included taxa (See Supplemental information). The morphological data matrix covers 36 taxa from all continents and spans the entire Cenozoic. The matrix contains 125 characters, 101 cranial and 24 postcranial, of which 96 are multistate; none was treated as ordered. The basal proboscidean Moeritherium was used as the outgroup.
We used Bayesian methods that allow for the simultaneous estimation of topology and divergence times, including fossil taxa as tips in a tip-dating approach. Several studies have shown that Bayesian methods combining molecular and morphological data and including fossil taxa, even when missing data is present in several taxa, are useful for obtaining accurate phylogenies (Wiens, 2003(Wiens, , 2009Wiens et al., 2010;Guillerme and Cooper, 2016;Mongiardino Koch et al., 2021). Specifically, we used a Birth-Death with Serial Skyline Sampling (BDSKY) tree prior (Stadler et al., 2013) with constant sampling over time, thus analogous to a BDSS tree prior (Stadler, 2010) or a fossilized birth-death (FBD) model (Heath et al., 2014). This analysis was performed with the software BEAST2 v2.3.2 (Bouckaert et al., 2014) using BEASTMasteR (Matzke, 2015) for the assembly of the xml file. For the morphological data, Lewis's (2001) Mkv model was used for the analysis. Variable gamma rates were favoured over equal rates (log BF > 200; stepping-stone sampling) after the implementation of the stepping stone method. Regarding the molecular data, due to limitations of the available models, the GTR nucleotide substitution model was used. A lognormal uncorrelated relaxed clock model was used for the analysis and uninformative uniform distribution [0,10] hyperpriors were used ll OPEN ACCESS iScience 25, 103559, January 21, 2022 iScience Article for the mean and the standard deviation. Additionally, birth and death rate priors were set to uniform distributions [0,100] and the sampling rate prior was also set to a uniform distribution [0,10]. Considering the extant taxa sampled in our dataset, Rho was fixed to one.
Biostratigraphic ranges of fossil taxa are provided in Table S3. All fossil taxa were treated as tips and calibrated with uniform distributions to account for stratigraphic uncertainty (Wood et al., 2013). Also, specific basal nodes were calibrated similar to the previous molecular-only analysis. The analysis was run for 100 million generations with sampling every 1000 generations. Convergence and effective sample size (ESS) were checked with Tracer v1.6. Finally, MCMC results were summarized in an MCC tree discarding the first 25% of trees as burn-in.
Diversification analysis. We performed a diversification analysis over the total evidence phylogeny to investigate any possible shifts in diversification during the history of the clade. For this purpose, we used the R package RPANDA (Morlon et al., 2016). RPANDA allows for the detection of diversification shifts in regions of a tree that have distinct branching patterns, using a model-free approach to efficiently summarize the shape of a phylogenetic tree by its spectral density. In particular, we used the function spectR to detect the optimal configuration of shifts based on the eigenvalues associated with the phylogeny and characteristics associated with the spectrum of eigenvalues. Once the optimal shift configuration was obtained, we used the BICompare function to map the shifts on the phylogeny, and chose the best configuration for the previously obtained number of shifts comparing the BIC values. Considering that our total evidence phylogeny was non-ultrametric due to the inclusion of fossil taxa, we used the argument ''non-ultrametric'' for the analysis. This approach allowed us to obtain the optimal number of shifts recognized for the data and also select the best position for those shifts based on a statistical approach.
Historical biogeography. We also performed a biogeography analysis using the currently available data to trace the evolutionary history of proboscideans geographically and chronologically. Reconstruction of ancestral distributions was based on the time-calibrated phylogeny. Because biogeographic models assume that processes happen at the species level, and using genera ranges at tips could result in invalid results, our analysis used the ranges of the most widely distributed species in each genus (Table S3). This approach represents clear limitations for the results of the analyses, mainly regarding the recognition of small-scale events and, especially, those occurring at the species level. In this regard, it should be made clear that our goal is focused on reconstructing the deeper and most important biogeographic events, therefore ignoring the potential events that could have occurred in relation to specific extinct species within genera (other than the most widely distributed one for the respective genus) or other yet-unknown extinct species. Therefore, our analyses focused on evaluating the most probable areas of origin of the different proboscidean clades, with special emphasis on their biogeographic distribution in the Americas. We defined six geographic areas: Africa, Europe, Asia, North America, Central America, and South America. The biogeographic distribution of proboscideans was modelled in the R package BioGeoBEARS (Matzke, 2013) using the DEC model, since this model permitted speciation via the biologically relevant mechanisms of widespread vicariance and subsequent sympatry (Givnish et al., 2016). Each taxon was assigned to a set of these areas. The maximum number of areas at a given node was set to five, which is the range of the taxon present in the largest number of areas (Gomphotherium). Models were time constrained at 22 Ma to account for the physical connection between Africa and Eurasia after 22 Ma. A constraint considering the passage of Proboscidea to South America was also modelled for two ages: 3.5 Ma for a late passage ( Figure S4) and 9 Ma for a possible early passage ( Figure S5).